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Abstract 



A two-dimensional Yukawa liquid is studied using two different nonequi- 
librium molecular dynamics simulation methods. Shear viscosity values in 
the limit of small shear rates are reported for a wide range of Coulomb cou- 
pling parameter and screening length. At high shear rates it is demonstrated 
that this liquid exhibits shear thinning, i.e., the viscosity n diminishes with 
increasing shear rate. It is expected that two-dimensional dusty plasmas will 
exhibit this effect. 
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Many-particle systems characterized by the Yukawa potential include a variety of physical 
systems, e.g. dusty plasmas, charged colloids, astrophysical objects and high energy density 
matter. The Yukawa potential 0(r) oc Q exp(— r/\n)/r models a Coulomb repulsion that is 
exponentially suppressed with a screening length Ad- Yukawa systems behave like a liquid 
when the temperature exceeds a melting point which depends on Q, X D , and particle spacing, 
e.g. [1,2]. 

Transport parameters of Yukawa systems - the diffusion coefficient [3], the shear viscosity 
[4-6] and the thermal conductivity [6,7] - have mainly been calculated for 3D systems, but 
there is now an increasing interest in 2D settings. For example, in dusty plasma experiments, 
charged microspheres suspended as a monolayer in a gas discharge make a 2D Yukawa 
system. By creating a shear flow in such a particle suspension, the viscosity was measured 
in recent experiments using 2D suspensions [8] and quasi-2D suspensions consisting of a few 
monolayers of charged microspheres [9]. The transport properties of such ultrathin liquids 
are also of interest as macroscopic analogs of molecular flow in nanoscience applications [10]. 

Transport coefficients are meaningful if they are part of a valid "constitutive relation" 
between the gradients of local variables and fluxes. For shear viscosity r/, the constitutive 
relation j y = —r][dv x (y)/dy] relates a momentum flux j y to the velocity gradient dv x (y)/dy, 
which is also termed the shear rate. In a non-Newtonian fluid, rj may vary with the velocity 
gradient, whereas in Newtonian fluids it does not. In particular, if r\ diminishes as shear is 
increased, the fluid is said to exhibit "shear thinning". This occurs in simple liquids [11], as 
well as in complex mixtures such as foams, micelles, slurries, pastes, gels, polymer solutions, 
and granular flows [12]. Recently, experimenters have claimed to observe shear thinning in 
dusty plasma liquids [9]. These reports motivate our simulations to search for the presence 
of shear thinning in 2D Yukawa liquids. 

Subsequent to the experimental measurement of viscosity in a 2D dusty plasma [8], a 
2D molecular dynamics simulation was used to obtain the shear viscosity from the Green- 
Kubo relations [13]. In this Letter we will go beyond the results of Ref. [13], which were 
performed for equilibrium conditions, by using here nonequilibrium simulations to search for 
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non-Newtonian behavior under conditions of a high shear rate. We will also compute the 
viscosity over a wider range of T and k. 

Our simulations use a rectangular cell with edge lengths L x and L y and periodic bound- 
ary conditions. The number of particles is between N = 990 and 7040. The system 
is characterized by dimensionless parameters T = Q 2 /47re aA; B T and k = a/ Ad, where 
a = (1 /nix) 1 / 2 is the Wigner-Seitz radius, with n being the areal density. Additional pa- 
rameters include the thermal velocity v = (2k B T '/m) 1 / 2 and the 2D analog of the plasma 
frequency uj p = (Q 2 /2ireoma 3 ) 1 / 2 , the shear rate 7 = dv x /dy, and its normalized value 
7 = (dv x /dy)(a/vo). Two types of molecular dynamics techniques are applied for the stud- 
ies of the shear viscosity. 

Method 1 reverses the cause-and-effect picture customarily used in nonequilibrium molec- 
ular dynamics: the effect, the momentum flux, is imposed, and the cause, the velocity gra- 
dient (shear rate) is measured in the simulation [15]. Momentum in the liquid is introduced 
in a pair of narrow slabs A and B, which are centered at y — L y /4 and 3L y /4, respectively. 
At regular time intervals r we identify the particles in slabs A and B having the highest \v x \ 
in the positive and negative directions, respectively. We then instantaneously exchange the 
v x velocity component of these two particles without moving the particles. This artificial 
transfer of momentum between slabs A and B (which is accomplished without changing the 
system energy) produces a velocity profile v x (y), the slope of which can be controlled by the 
frequency of the momentum exchange steps. The equations of motion 



where r = (x, y), p = (p x , p y ) are the positions and the momenta of particles, m is their mass, 
and Fj is the force acting on particle i, are integrated by the velocity Verlet algorithm. 

Method 2 simulates a planar Couette flow, which is established by the Lees-Edwards 
periodic boundary conditions which result in a homogeneous streaming flow field in the 
simulation box: (v x ) = 7(2/ — L y /2), where () denotes a time average. The system is 
described by the Gaussian thermostated SLLOD equations of motion [16]: 



dTj _ pi 

dt m 



dpi = 
dt 



(1) 
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where p = (p x ,P y ) is the peculiar momentum of particles, x is the unit vector in the x 
direction, and a is the Gaussian thermostatting. multiplier. The above set of equations is 
solved using an operator splitting technique [17]. 

In contrast to method 1, method 2 results in a homogeneous shear field and a constant 
temperature within the whole simulation box. Thus arbitrarily high shear rates may be 
established without the need of considering any effects of temperature gradients on the 
viscosity. 

In both methods the pairwise Yukawa interparticle forces are summed over a /t-dependent 
cutoff radius, using the chaining mesh technique. (The force due to particles at the cutoff 
radius is ~ 1CT 5 smaller compared to that due to the nearest neighbors.) Both methods 
neglect any neutral gas drag, which has been observed to alter the velocity profiles in ex- 
periments [8], i.e. they model an atomic system where momentum transfer is dominated by 
Coulomb collisions [1]. 

Method 1 has the advantage that is resembles more closely the experimental conditions, 
although the procedure for applying shear in the simulation involves no introduction or re- 
moval of energy from the system. In the experiment [8] shear is applied via an external 
introduction of both momentum and energy in a boundary slab while energy is simultane- 
ously removed elsewhere by frictional dissipation. Method 2 represents a well-established 
technique for measurement of viscosity at arbitrary steady, as well as temporally varying 
shear rates. Although it has little connection to the conditions found in the experiment, 
it has been demonstrated to be an efficient technique to investigate shear thinning [11,16]. 
Thus we apply this method for the studies of this latter effect. 

Near equilibrium (small 7) shear viscosity values have been obtain using both techniques. 
In method 1 this is done at the lowest practical shear rate, where dv x /dy is uniform between 
slabs A and B. We calculate r/ eq from 



j y \ = i] cq dv x (y)/dy = Ap/2t shn L y , 



(3) 



where Ap is the total x-directional momentum exchanged between slabs A and B during 
the simulation time £ sim [15]. In method 2 the off-diagonal element of the pressure tensor is 
measured during the course of the simulation: 



In method 1, the spatial profiles for temperature and velocity, Fig. 1, develop self- 
consistently in response to the perturbation applied by introducing momentum in slabs 
A and B. We use method 1 only for small perturbations, so that the velocity profile 
has a linear gradient and the temperature is isotropic, with T x = T y , where T x>y = 
(m/Njk-B)Y^i=!:i([vix,y(t) — Vj^j] 2 }- The index i runs over the Nj particles in slab j. We 
verified that v]y~ is negligibly small. 

Obtaining reliable results for r\ at small 7 requires a simulation duration of typically 
uOpt ~ 10 4 - 10 5 for both methods. The required time step is smallest, and the simulations 
are most costly, at low Y. In method 1, system size effects are expected to appear when 

(i) particles traverse the simulation box without significant interaction with the others, or 

(ii) the compressional sound wave transits the box in a shorter time than the decay time t c 
of the velocity autocorrelation function. For (i), for our most demanding condition (small 
size N = 990 and high temperature T = 1) a particle moving at the thermal velocity would 
transit the cell in a time oo p t « 57 if it were undeflected by collisions. We find that the 
decay time is short enough, u p t c ~ 5 - 10, even for the smallest T values of interest. Thus 
we expect no "ballistic" trajectories across the entire simulation box. For (ii), the sound 
speed [14] at k — 1 is v — dw/dk ~ au p , and the wave's transit time At across a box with 
length L y is uj p At ~ L y /a = 57 for our most demanding case, N = 990 particles. Thus, we 
find both criteria fulfilled for a "sufficiently large" system. Method 2 is known to produce 




(4) 



where = — rj = (xij,yij), and the shear viscosity is obtained as 



r] = \mi(P xy (t))h. 

t— too 



(5) 



5 



accurate results even for small number of particles simulated [16]. We verified that the 
results obtained from both methods did not depend significantly on N. 

Figures 2(a) and (b) illustrate particle trajectories in simulations based on method 2, 
for conditions r = 10, k — 1 at a shear rate 7 = 0.2, and for T = 100, k = 1, 7 = 0.05, 
respectively. 

Our results for r] cq as a function of T, for different values of k are plotted in Fig. 3(a). We 
find a good agreement with the earlier equilibrium MD simulation of Ref. [13]. In contrast 
with most simple liquids, which have a viscosity that varies monotonically with temperature, 
a prominent feature of the viscosity of the present system is a minimum (e.g. at T = 20 
for k — 1) , which has been noted previously in both OCP (one-component plasma) and 
Yukawa liquids. The shape of the i] eq (T) curve can be explained by the prevailing kinetic 
and potential contributions to the viscosity at low and high values of T, respectively. The 
near-equilibrium shear viscosity values obtained with method 2 for k = 1 are also displayed 
Fig. 3(a). We find an excellent agreement between the results of methods 1 and 2. 

Similar to what was observed in [5] for 3D Yukawa liquids, we find that the near- 
equilibrium viscosity i] cq obeys a scaling law as demonstrated in Fig. 3(b), where viscos- 
ity has been normalized by r/ E = mnui-^a 2 . The Einstein frequency ue depends on k, and 
we computed it from Eq. (7) of Ref. [14] using pair-correlation functions measured from 
our simulations. The horizontal axis is a normalized temperature T' = T y /T m = T m /T, 
where T m and T m are melting-point values reported in Ref. [2]. Using these normalizations, 
the data fall on the same curve, demonstrating the existence of a scaling law for the 0.5 
< k < 2.0 range of the screening parameter. We note that for this purpose we found u E 
was more significant than u p . The near-equilibrium viscosity is fit by an empirical form 
VcJve = aT + b/T' + c with coefficients: a = 0.0093, b = 0.78 and c = 0.098. 

A shear-thinning effect is revealed in Fig. 4(a), which shows that rj diminishes significantly 
as the shear rate 7 is increased. In other two-dimensional systems the reduction in r], as 
compared to the value at small shear, was observed to vary as the square root of 7 [11]. 
We find that this scaling also occurs for the Yukawa system, as indicated by data that fall 
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nearly on a straight line in Fig. 4(a) for 7 > 0.2. At smaller shear rates, 7 < 0.2, however, 
the shear thinning effect is less profound and the liquid is more nearly Newtonian, especially 
for large T. Results are shown for 7 > 0.01, which we found to be reliable, whereas at lower 
7 method 2 yielded noisy data even for very long simulations. 

Because viscosity arises from both kinetic and potential contributions, Eq.(4), we eval- 
uate which of these contributions is most responsible for the observed shear-thinning effect 
in Fig. 4(b). Recall that for equilibrium conditions, the kinetic term dominates for T <C 
20 and the potential term dominates for T S> 20. Here, we find that for non-equilibrium 
conditions, as the shear rate increases the reduction in viscosity is mostly due to a reduction 
of the kinetic contribution at low T and a reduction of the potential contribution at high 
T. In other words, at extreme values of T, it is the same term dominating the equilibrium 
viscosity that also dominates the shear thinning effect. At an intermediate value of T = 20, 
however, where the equilibrium viscosity has a minimum and the two equilibrium contribu- 
tions are comparable (for k — 1), we find that it is mostly a reduction of the kinetic term 
that accounts for the observed shear thinning effect. 

In summary, we have calculated the shear viscosity coefficient of 2D Yukawa liquids in 
a wide domain of parameters T and k using two different molecular dynamics approaches. 
The small shear rate calculations confirmed that the two techniques used are consistent 
and yielded r\ values in fair agreement with equilibrium MD calculations [13]. The small 
shear rate data were found to obey a universal scaling: r\ normalized by the Einstein fre- 
quency was found to depend only on the reduced temperature (ratio of the temperature to 
melting temperature). The high shear rate simulations based on method 2 unambiguously 
demonstrated a non-Newtonian behavior of the Yukawa liquid: r\ was significantly reduced 
for these conditions in a manifestation of shear thinning, except at the lowest shear rates 
where the liquid is more nearly Newtonian. Regimes of the plasma coupling parameter were 
identified to distinguish whether the kinetic or potential contribution to the shear viscosity 
is primarily responsible for the shear thinning effect. 
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FIGURES 

FIG. 1. (a) Velocity profiles v x (y) obtained from method 1 for different frequencies (1/r) of 
momentum exchange steps, and (b) T y (y) temperature profiles for the same conditions. V = 100, 

K = 1. 

FIG. 2. (a) Trajectories of particles in the simulation based on method 2. (a) r = 10, n = 1, 
at a shear rate 7 = 0.2, time of recording : w p AT = 5.0; (b) T = 100, k = 1, 7 = 0.05, w p AT = 
23.6. The shear field is = 7(2/ — L y /2), i.e. there is no flow at y = L y /2. N = 1020. 

FIG. 3. (a) Shear viscosity at near-equilibrium conditions, rj eq , obtained from method 1, at the 
simulation's lowest practical shear rate, and normalized by rjo = mnuj p a 2 . For comparison, data are 
shown from the Iowa equilibrium MD simulation and from simulations based on method 2 (SLLOD), 
in the limit of small shear rates, at k = 1. N is the number of simulation particles, (b) A scaling 
law is demonstrated by normalizing the data in (a) using tje = mnujpa 2 and T' = T y /T m where 
T m is the melting temperature. The thick line is an empirical fit of form ^ eq /??E = oT' + b/T' + c. 

FIG. 4. (a) Shear viscosity as a function of normalized shear rate 7 for selected values of the 
coupling parameter V. (b) Potential (filled symbols) and kinetic (open symbols) contributions to 
the shear viscosity, k = 1. 
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This figure "fig2.gif.gif" is available in "gif" format from: 



http://arXiv.org/ps/cond-mat/0603667vl 
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